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PREFACE 


The objective of this program was to demonstrate the capability of the 
Westinghouse image registration techniques to perform a geometric quality 
assessment of muLtispectral and multitemporal image pairs. Registration 
measurements were performed by computer processing of images supplied 
by NASA-Goddard under Contract NAS5-20947. The program was conducted 
by the Westinghouse Systems Development Division, Baltimore, under the 
cognizance of Mr. Bernard Peavey of NASA-Goddard, The Westinghouse 
program manager was Dr. Glenn E, Tisdale, Computer prograrnming was 
performed by Mr, David T. Bissell, 

It was concluded that the Westinghouse techniques are capable of provid- 
ing multispectral or multitemporal registration of LANDSAT imagery to 
accuracies of a small fraction of a pixel. They are insensitive to the choice 
of registration area, provided terrain detail is present. They may be im- 
plemented at megapixel-per- second rates using commercial minicomputers 
in combination with a special-purpose digital preprocessor. 

It is recommended that implementation of the techniques as a quality 
control tool be initiated by conducting a brief systems analysis. 
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1. INTRODUCTION 


The prodigious data gathering capability of the earth resources satellites 
offers wide application to many fields. It also provides the challenge to 
extract useful information with a minimum of data processing. One solution 
to this problem is to follow the initial assessment of a particular area with 
the further detection of changes only. However, change detection requires 
detailed registration between the original (reference) image and the more 
recent sample. Because of differences in spacecraft geometry and condition, 
these multitemporal images, although visually similar, will not precisely 
register. Unless the registration can be held to within a small fraction of 
an image element (pixel), the change detection process will generally 
produce errors or false alarms. 

A solution to the registration problem lies In the modification of the 
geometry of one image so as to correspond to the other. This could be 
accomplished by visual comparison between corresponding well-defined 
image points or objects. Aside from the effort involved, this process 
suffers from the difficulty of defining corresponding points to within a small 
fraction of a pixel. A more satisfactory answer lies in the comparison of 
corresponding image areas using automatic correlation techniques. Im- 
proved accuracy is possible because the displacement values depend on an 
array of image points. Random errors inherent in the location of a single 
point are averaged out by the correlation process. 

Conventional correlation techniques operate by comparing sample image 
areas on a point-by-point basis. This process has the following disadvan- 
tages: 


a. Because all of the image points are involved on an equal basis, 
the process is inefficient from an information handling point of view. 

Most image points (in areas of uniform density, for example) carry little 
or no information. As a result, processing systems require large band- 
width and storage capability. 

b. Coi'relation of image density values is degraded by density varia- 
tions caused by ambient illumination, seasonal changes, and differences in 
sensor frequency (multi spectral correlation). 

c. Initial acquisition between images is inefficient and slow if mis- 
registration is sizeable, or if differences exist in relative orientation and 
scale. 

d. Individual correlation measurements provide no indication of the 
direction in which a match condition will be located. 

As an outgrowth of a research program in pattern recognition, Westing- 
house has developed digital image registration techniques that resolve the 
problems stated above. These techniques were applied originally to auto- 
matic target recognition tasks, where very significant differences between 
images must be accommodated. As a result, they operate on the global 
pattern of an image area, as opposed to the specific density samples. 

The objective of the reported program was to investigate the feasibility 
of using the Westinghouse image registration techniques as a quality control 
tool, offering rapid assessment of the geometric comparison between image 
pairs. The program was initiated in 1975 by Mr. Bernard Peavey of 
NASA-Goddard, following the description of the techniques by Westinghouse. 
The goals were to demonstrate and measure registration performance on 
IA.NDSAT images supplied by NASA and to establish the feasibility of im- 
plementing the capability at megapixel -per- second rates, • 
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The report will proceed in Section 2 with a detailed review of Westing- 
house image registration techniques. The test program to measure 
registration performance is described in Section 3, Implementation at 
me gapixel-per- second rates is discussed in Section 4, Finally, conclusions 
and recommendations are made in Section 5. 
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2. WESTINGHOUSE TECHNIQUES FOR IMAGE REGISTRATION 


A fully automatic image registration technique should provide for initial 
acquisition of the match between the images as well as the precise com- 
parison of their individual points. The initial acquisition should be possible 
even though tlie images differ in relative orientation, position, and scale. 
Some existing correlation techniques require close manual alignment be- 
tween the images before detailed comparison is possible. The Westinghouse 
techniques, based on the recognition of global image patterns, provide auto- 
matic acquisition capability, as well as precise final alignment of image 
points. They are based on a company- sponsored development started in 
1965. U. S, patent No. 3, 748, 644, entitled '^Automatic Registration of 
Points in Two Separate Images, ” was issued to Westinghouse on July 24, 
1973, to cover this approach. In addition to the concepts covered in the 
patent, Westinghouse has developed proprietary methods for digital image 
processing that are incorporated in the registration operations. 

Since the digital image preprocessing methods are basic to the imple- 
mentation of the Westinghouse registration technique, they are described 
in paragraph 2. 1. In particular, this paragraph describes a method for 
image preprocessing that greatly reduces image data bandwidth and extracts 
key image information in the form of edge contours and overall statistics. 
Next, paragraph 2. 2 is a description of how these preprocessor outputs are 
used in the general case of initial acquisition between images when limited 
prior knowledge is available with respect to their relative orientation, 
position, and scale. After completion of initial acquisition, or if sufficient 
prior information is available to permit initial positioning, precision regis- 
tration or tracking may be accomplished at high speed, using a technique 
described in paragraph 2. 3. 
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2, 1 Digital Image Preprocessor 

Image registration operations are accomplished by a hybrid digital image 
processor, consisting of a special-purpose preprocessor (to extract key 
image data) followed by a programmable general-purpose processor (to 
extract registration vectors), A description of processor hardware is 
contained in Section 4, In this paragraph, the operations in the special- 
purpose preprocessor are described. Registration operations are described 
in paragraphs 2. 2 and 3. 

The function of the preprocessor is to extract the features required for 
registration computations from the gray level image samples. These 
features are straight-line contours of density gradient. In effect, the pre- 
processor produces a line drawing of the image. At the same time, the 
data bandwidth is reduced by one or two orders of magnitude, depending on 
the selected thresholds for minimum gradient detection and line length. In 
addition, the preprocessor routinely extracts statistics associated with a 
uniform matrix of image areas. During this program, these uO -ituxal 
statistics were used in the selection of appropriate prefilters, as dis- 
cussed in paragraph 3, 3, 

Operation of the preprocessor is on a line-by-line basis with respect to 
the input image. Therefore, video data may be handled directly. Further- 
more, storage requirements in the preprocessor are limited to single lines 
of data only. The preprocessor comprises a series of shift registers. 

Data flow through the line extraction portion of the preprocessor is shown 
in Figure 2-1. 

We will demonstrate the individual preprocessor operations with the use 
of LANDSAT images 1703-17590 (Reel 2) and 1739-17575 (Reel 2) repre- 
senting location N 36 deg, W 120 deg (near Fresno, California) on June 26, 
1974, and August 1, 1974, respectively. Portions of these images (MSS 
Band 4) are displayed in Figure 2-2. In particular, we will concern our- 
selves with the windows A and A*, measuring 125 x 125 pixels in size. 
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Figure 2-1. Flow of Data Throxigh the Preprocessor 

Gray level density values for these two areas are shown by Figure 2-3. 

The numbers represent the density to the nearest of 16 levels. Numbers 
above 9 are indicated with an overprinted "/”. Large numbers represent 
dark areas. 

2.1.1 Image Prefiltering 

Before preprocessing the images for line extraction! spatial filtering 
operations may be performed on the gray level data for purposes of noise 
smoothing or equalization of resolution in both dimensions. It was found 
desirable to provide filtering between the sensor scan lines of the LANDSAT 
data, as explained in paragraph 3. 3. Filtering between lines requires 
storage of one line of data samples in addition to those shown in figure 
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Figure 2-2. Portions of LANDSAT Images - 
of Fresno, California for June 26 
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2, 1. 2 Gradient Extraction 

To generate the straight~line contours (subsets) of the images it is first 
necessary to compute the two-dimensional gradient (derivative) at each 
image point. As shown in Figure 2-4, this is done with a four-pixel window, 
which scans across the image in a raster format. The gradient amplitude 
and direction are approximated as shown. The gradient direction (0 to 360 
degrees) is quantized to one of 16 discrete directions, as depicted in the 
diagram. To suppress the areas of negligible gradient activity (containing 
no significant contour or edge information), a threshold is applied to the 
gradient amplitude. Gradients with amplitudes less than the threshold are 
set at 0. A display of gradient directions for windows A and a' in Figure 
2-2 is shown in Figure 2-5. The numbers indicate direction from light to 
dark in accordance with Figure 2-4, Numbers above 9 are overprinted 
with a "/". 

2, 1. 3 Gradient Maximizing 

After gradient thresholding, the edges are generally still too wide for 
subset generation. Therefore, a gradient thinning operation is performed 
that basically “skeletonizes*' wide edges. 

The algorithm uses a raster scanning window containing a gradient cell, 

X, and four of its nearest neighbors. The neighbors with colinear gradient 
directions are compared to X, The largest gradient amplitude is then re- 
tained, This procedure is repeated sequentially for each gradient point in 
the image. The maximized gradient display for images A and a‘ is shown 
by Figure 2-6; compare this with Figure 2-5, 

2,1.4 Subset Generation 

Subset generation is accomplished by "growing" a line formed by adjacent 
parallel gradients. As before, a 5- cell scanning window is employed (see 
Figure 2-7). The new data point is labeled cell X. Its four neighbors are 
examined to find those containing a parallel (within a tolerance) gradient 
direction. If one is found, then the neighbor is dropped as a line endpoint. 
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GRADIENT 

COMPONENTS 

A X = B -C 
A Y = D ■ A 


GRADIENT AMPLITUDE = MAX.( Ax, A Y) + ’/j MIN. ( A X, A Y) 

1 Ay 

GRADIENT DIRECTION = [ TAN"^ ), QUANTIZED INTO 1 OF 16 DIRECTIONS FOR 0-> 2 
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Figure 2-6. Maximized Gradients for the Windows in Figure 2-3 
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and X becomes the new endpoint. Neighbors that are colinear to the gradient 
of X are excluded to prevent false lines from forming. The operational 
cycle of the subset generator is diagrammed in figure 2-7, 

The subsets derived from images A and A* are shown in the graphical 
plots in Figure 2-S. The subsets are numbered consecutively according to 
their lowest position on the plots. The polarity (location of light and dark 
sides) for each subset is coded by the presence or absence of stars at both 
ends. In general, the absence of the star indicates a change from light to 
dark in the clockwise direction, for an axis of rotation about the uppermost 
endpoint. For horizontal subsets, the a^s of rotation is placed at the right 
end of the subset, 

2,1.5 Texture Statistics 

A second preprocessor function is the collection of statistics over a 
matrix of square image areas. The measurements amount to a textural 
analysis of the windows. They include histograms of gray density, gradient 
amplitude and direction, and subset length. Average values for gray 
density and gradient amplitude are also computed. 

2, 2 Initial Acquisition 

The 'Westinghouse image processing techniques can achieve registration 
between two line-by-line scanned images of the same area totally independent 
of relative orientation and independent of position and scale over reasonable 
ranges. The registration output includes a measure of image similarity 
between the images, their relative scale and orientation, and the transfor- 
mation of points between images. These steps are carried out in the sequence 
shown in Figure 2-9. 

With reference to Figure 2-9, video signals from both images are digitized 
and preprocessed to extract image density contours or subsets, as described 
in paragraph 2. 1, The subsets are defined by their endpoints. Features are 
formed from the subsets that consist of the geometric relationships between 
pairs of subsets and provide properties that are invariant, regardless of 
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Figure 2-8, Line Contours (Subsets) for the Windows in Figure 2-3 



Figure 2-9. Registration Steps for Initial Acquisition 

relative position, orientation, or scale. For example, these properties can 
include the relative angles between subsets. These features are compared 
between images, with the initial result being a measure of similarit 7 asso- 
ciated with the registration, as well as a measure of their relative orienta- 
tion and scale. When the similarity match has been accomplished, endpoints 
on the corresponding contours in the images can be related. In particular, 
the center point of one image can be located on the second image. 

The computations associated with a pair of subsets are listed in Figure 
2-10, Subsets AB and CD in image A correspond to A*B' and C’D' in image 
B. A geometric feature is formed between the pairs by connecting endpoints 
A and C and A' and C. The angles between the subsets and the connecting 
lines, 7-, 7_, and t' and 7i will be equal if figures ABDC and A’B'D*C* are 
similar. Furthermore, their values are invariant with translation, scale, 
and orientation of these figures in the plane. Accordingly, these angles are 

















Final Transforms 


X‘ = SVS [(X + Xq) cos - {Y-S-Yq)S!N A■^•] 

Y' = S7S OX + Xq) SIN A<^ + iY+Yg) COS 

74-0074-VA-336 

Figure 2-10, Computations in the General Case of Initial Acquisition 
Ijetween Images with Uncertainty of Position, Orientation, and Scale 




the first elements in the geometric features used to achieve automatic 
acquisition^ Along with these angles* the features also include the orienta- 
tion, Ip, <f>^, of the connecting lines relative to a reference axis in the 
image, and the lengths, S and S’, of the connecting lines. They also include 
the subset endpoints. 

Initial acquisition is obtained by matching pairs of features between the 
images with respect to 7, and by clustering the results with respect to rela- 
tive orientation and relative scale. If the two images are identical except 
for translation, rotation, and scale, all the features will cluster at a point. 

If they have relative distortion, the cluster will be dispersed. The amount 
of acceptable distortion is controlled by an adjustable threshold. 

The center of the cluster provides values for the relative orientation and 
scale of the two images. These values may be used to determine a transfor- 
mation between the x and y coordinates for the two images. Since corres- 
ponding features contain endpoint coordinates, each pair can be used to 

calculate the nominal coordinates of translation, x , y , between the 

o o 

images. These are again clustered to obtain an average over the images, 

and to eliminate anomalies caused by symmetrical figures such as rectangles. 

The cluster center in the x , y plane is used with the S values to com- 

o o 

plete the final transformation from any point in image B to the corresponding 
point in image A, 

Where correction for image warp is desired {for example, in change de- 
tection operations), registration is carried out on a matrix of windows having 
stafficient density to accommodate the warp function. 

If one or more parameters, such as relative scale, are known, then the 
process can be carried out with less computation, and with improved pre- 
cision, than would otherwise be the case. 
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2. 3 Precision Registratiojo . 

The process of precision registration can begin when the initial acquisition 
between images has been completed, or it can be carried out directly if 
relative scale, orientation, and position between the images results in 
differences between them amounting to only a few pixels. Precision regis- 
tration makes use of preprocessor subsets, as did the initial acquisition 
process. However, the computations are greatly simplified. In fact, a 
positional search, with scale and orientation known, is generally most 
readily accomplished by performing a matrix of precision registration 
measurements. These operations are demonstrated below. 

The precision registration process involves the matching of corresponding 
subset endpoints between the two image areas and the computation of their 
average relative displacements by summation over the areas. With reference 
to Figure 2-11, windows are erected in one image around the location of 
subset endpoints from the other image. If subsets of appropriate direction 



Figure 2-11, Matching of Subset Endpoints in the Precision 

Registration Process 
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terminabe within the windows, the displacement of the endpoint relative to 
the window center is computed. For the window at upper left in Figure 
2-11 (Point 1), the displacement is Ax = 0, Ay = 0. For the window at 
lower right (Point 2), the displacement is Ajc = +1, Ay = +1. 

Over the entire image area, the coordinates for the average displacement 


will be 


Ax = 


Ay = 


2 Ax 


( 2 - 1 ) 


where ti is the number of matching endpoints. 

The registration algorithm includes a procedure fOi handling the occurrence 
in the window of more than one stibset. In the event that no matching subset 
is found, there is no contribution to the summation. 

The error in the displacement of computation of equation 2-1 depends on 
the positional uncertainty of the subset endpoints relative to the actual 
ground characteristics they represent and the number of measurements, n, 
over the area, as follows; 


Error = 


Positional Uncertainty 


( 2 - 2 ) 


The positional uncertainty includes the effects of noise in the sensor system, 
the quantizing effects of positioning to the nearest pixel, and the character- 
istics of the preprocessor in assigning endpoints to the subsets. The summa- 
tion of these effects is expected to produce a positional uncertainty of about 
2 pixels with the LANDSAT imagery, and in a random direction. 

The positional uncertainty in equation 2-2 involves a pair of endpoints, 
one in each image. Accordingly, it is larger than the estimate for one end- 
point as stated above. If we designate the positional uncertainties for the 
endpoints as a and b, and the angle between them as 0, then the net indicated 


uncertainty, u, will be given by 
2 2 2 

u = a +b 2ab cos 0, 


(2-3) 
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from simple geometry. Since the angle between them is assumed to take 
any value at random, the average value for the term on the right is zero. 
Accordingly, a good estimate for the combined positional uncertainty for a 
pair of endpoints is the rms value for two endpoints, or about 2. 8 pixels. 

If 100 endpoint pairs are involved in a displacement measurement, in 
accordance with equation 2-2, the expected error will be 0. 28 pixel. If 
400 pairs are available, the error would drop to 0, 14 pixel. If 800 were 
available, the expected error would be 0. 10 pixel. An increase in image 
area usually results in a corresponding increase in the number of available 
endpoints. Therefore, if an area of 125 pixels on each side produces a 
relative displacement accuracy of 0. 28 pixel, an area of 350 pixels on a 
side should produce an accuracy of 0, 10 pixel. 

In Section 3, accuracy measurements are discussed that can be compared 
with the above estimates. 

2. 3. 1 Computation of Relative Displacement 

The characteristics of the displacement computation are readily 
visualized if a matrix of registration measurements is made in the vicinity 
of the position of best match. Such a matrix is useful for initial acquisition 
if relative translation is present, but orientation and scale differences are 
minor. However, few of these computations are ordinarily required when 
acquisition has been made. 

A matrix of measurements for the windows A and A' of Figure 2-3 is 
described in Figures 2-12, 2-14, and 2-15. The matrix includes 20 relative 
positions in both x and y, or a total of 400 computations. The matrix of 
Figure 2-12 shows the number of matching endpoint pairs (divided by 100). 

In this computer printout, the numbers along the axes relate to the coordi- 
nates of the windows on the NASA digital tapes. The crossed lines locate 
the position of best match, which will be subsequently computed. The box of 
six pixels on a side corresponds to the size of the match window (see 
Figure 2-11). 
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Figure 2-12. Matrix of the Niamber of Matching Subset Endpoints 
(-riOO) for Windows A and A’ of Figure 2-3 

Note that in the vicinity of the best match position, the number of match- 
ing pairs rises well above surrounding values. The increased values indicate 
an increased correlation for these locations. Values for best match repre- 
sent approximately one- third the available subset endpoints. 

Some matches will generally occur for any relative position between the 
images. The actual number will depend on the si:se of the match window, 
the density of subset endpoints, and the specific algorithm for match 
acceptance. In general, a compromise will be required to balance an ade- 
quate level of matches to provide for accuracy, with a limited number of 
random pairings. 

To appreciate the method for computation of the match position, it is 
useful to visualize the effect on the offset computations of sliding one set of 
endpoints relative to the other. In Figure 2-13, this effect is shown for one 
axis. A match window of 6 pixels in width is assumed (seven match positions). 
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Figure 2-13. Indicated Displacement Between Images as a Function of 
Relative Displacement for a Single Image (Curve A) and Separate 
Images (Curve B) of the Same Area 


The effect of sliding one image against itself is shown as Curve A. For this 
curve, random matches will occur until the relative displacement is equiva- 
lent to -3 pixels, at which point a match would be found between all end- 
points but with an indicated displacement of -3 pixels. For the next 6 pixel 
positions, a linear curve of displacement would be described, as shown. 
When the relative displacement becomes +4, however, the output would again 
be the result of random matches. 

For different images of the same area, a result similar to Curve B would 
be obtained. The values for indicated displacement are reduced by the 
positional uncertainty of the endpoints in each image. At the edges of the 
window, for example, this uncertainty will result in the loss of some match- 
ing pairs, while indicated 'displacement for others will appear at low values. 


I 
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Calculation, of the relative displacement under these conditions is carried 
out by a determination of the point of crossing the relative displacement 
axis. This is obtained by linear interpolation between the nearest integral 
positions. 

Returning to images A and A', Figure E-14 provides a matrix of values 
for the indicated displacement in the x direction {along the sensor scan lines) 
The corresponding displacement in the y direction is given by Figure 2-15. 

In Figure 2-l6, the indicated displacement versus relative displacement 
from Figure 2-14 are plotted for both the vertical positions of 4 and 5. The 
similarity is obvious. Computation of the crossover point by linear inter- 
polation would yield: 


A 0.3Q 

^4 ^ 0. 30 + 0. 15 


= 0.67 


A 0. 36 

^5 ~ 0. 36 -i- 0. 16 " * ^ 


(2-4) 
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Figure 2-15. Matrix for Indicated Displacement in tne y Direction 
for Windows A and A' in Figure 2-3 

An average for these two values is 0. 68. Since we are interested in per- 
forming registrations at high speed, computation of displacement was 
carried out for the present tests by using for both x an.d y the values for 
indicated displacement which are increased in x and y by one unit. Under 
these conditions, the computation for Ax would become 


Q. 30 

0. 30 + 0, 16 


= 0.65 


(2-5] 


This value differs from the average from equation 2-4 by 0. 03 pixel. Since 
the number of matches (140) would indicate a positional accuracy of 0. 24 
pixel, the difference seems acceptable. 
















In the y direction the £ig*'r"z nre extracted from Figure 2-15 in the down- 
ward direction for horizontal positions -1 and 0. We have 


0.35 

“ 0. 35 + 0. 35 
0. 48 

^^0 ^ 0. 48 + 0. 29 


0. 50 
0. 63 


+ Ay^ = 0, 57 


(2-6) 


0.35 

^-1, 0“ 0. 35 + 0. 29 


s 0. 55 


The difference is 0,02 pixel in the y direction. 

In the test program described in Section 3, displacements were computed 
by linear interpolation across diagonal positions. 
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3. TEST PROGRAM 


The objective of the test program conducted for this study was to estab- 
lish the feasibility of the Westinghouse image registration techniques to 
register LANDSAT imagery, both multispectral and multitemporal, and to 
estimate the accuracy with which this could be accomplished. The following 
tests were completed; 

a. Multispectral registration between all combinations of the four 
MSS bands 

b. Determination of an appropriate prefilter for the preprocessor 
and selection of preprocessor thresholds 

c. Measurement of a distortion matrix between bulk multitemporal 
images 

d. Measurement of a distortion matrix between previously processed 
and corrected multitemporal images 

e. Quality measurements on specially processed multitemporal 
samples. 

Paragraph 3, 1 is a desci'iption of the methods for preparing and formatt- 
ing NASA digital tapes. The tests listed above are described in paragraphs 
3, 2 to 3, 6. An estimate of registration accuracy is computed in paragraph 
3. 7 from the data collected in the individual tests, 

3, 1 Preparation of Test Data 

All image data used in performing registration measurements were obtain- 
ed from NASA-Goddard on computer-compatible magnetic tapes, ^ They were 
processed at Westinghouse on a Uni vac 1110 system. 

*Tape formats are described inNASA-GSFC Report X-563-73-206, ‘'Gener- 
ation and physical Characteristics of the ERTS MSS System Corrected Com- 
puter Compatible Tapes, " by Valerie L. Thomas, July 1973, 


The first step in data preparation is to copy the tapes from 9-track 
800-bpi format to the 1110 system standard of 9-track 1600-hpi* Copying 
the tapes provides a dump of the header identification, and annotation 
records, hut it does not separate the four image bands. 

The next step is to extract windows of 125 by 125 pixels each at desired 
locations and to separate the four spectral bands. The choice of 125 by 125 
pixels provides sufficient image area to verify registration measurements, 
but with low individual processing times. However, other sizes may be 
selected if desired. Where a uniform matrix of windows is desired, these 
may be extracted automatically by the extra/ tion program. The output of 
the program is a magnetic tape containing a sequence of four windows for 
each of the four spectral bands at the desired locations. The density values 
for each image sample are complemented and reduced to six bits of resolu- 
tion for input to the registration logic. 

3. 2 Multispectral Registration 

In testing a registration method, the measured performance must be 
compared with some reliable standard, or ground truth. The ground truth 
may be obtained from 

a. An exact knowledge of sensor geometry 

b. An alternate, more reliable means for performing the measurements 

c. By averaging a series of unrelated measurements at the same points. 
The positioning accuracy and linearity of the LANDSAT images preclude the 
reliance on sensor geometry for images taken at different times (multi- 
temporal). However, the images from the four MSS spectral bands for the 
same location are presumably in precise registration. Accordingly, initial 
testing consisted of the computation of the six multispectral combinations of 
these bands. The results are displayed in Table 3-1. Three scenes, each 
125 by 125 pixels in size, were chosen from ERTS-E- 1757- 17574, located 

on the California coast at longitude 120 deg W, and latitude N 34 deg 30 inin. 
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3-1 


FOR MULTISPECTRAL IMAGES 


Relatives 

Displacement 

Components 


X 


Value of 
Displacement 

Matches 

Subsets 

0.30 

90 

60/190 

0.44 

71 

60/221 

0.66 

26 

60/49 

0.10 

159 

190/221 

0.69 

50 

190/49 

0.32 

85 

221/49 

0.24 

234 

168/274 

0.29 

185 

168/325 

0.17 

72 

168/184 

0.21 

243 

274/325 

0.76 

89 

274/184 

0.15 

243 

325/184 

0.04 

246 

166/196 

G.18 

57 

166/175 

- 

19 

166/143 

0.39 

65 

196/175 

- 

23 

196/143 

0.25 

161 

175/143 


- 0.02 
0.12 
Wo Lock 

j ^0.13 

No Lock 

■ 0.01 
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Scene 1 is a mountainous area, while scenes Z and 3 are urban in nature. 

In all tests, the x coordinate is taken along the MSS scan (east-west), while 
the y coordinate is between scans (north-south). 

Displacement components in Table 3-1 were calculated in the manner 
described in paragraph 2-3, For these multispectral samples, the values at 
the X = 0, y - 0 position in the displacement matrix should be 0, Actual 
values are tabulated in the third column of the table. A crossover point is 
obtained by selecting the adjacent diagonal indicated displacement value, 
which exhibits a change in the sign of both components, and by performing a 
linear interpolation. The resulting relative displacement components are 
shown in the fifth column. The amplitude of this displacement is also 
tabulated, together with the number of endpoint matches, and the number of 
subsets in the images from which the matches ai-.; obtained. 

It is of interest to plot the error as a function of the number of matches 
in the window. The displacement values of Table 3-1 are located in Figure 
3-1 in this manner, with the pair of spectral bands indicated beside each point. 






Also shown are curves that describe the relationship of equation 2-2 for 
values of positional uncertainty for a pair of endpoints of 2, 4, and 6 pixels. 
Except for image pairs including band 4, all measurements are within the 
curve labeled 4, If the assumption is made that the displacement components 
in Table 3-1 are samples from independent Gaussian distributions and that 
equation 2-2 holds, then a maximum likelihood estimate of relative positional 
uncertainty is 3. 43 pixels, using all the data. This implies 3, 43/s/2 = 2. 43 
pixels positional uncertainty for a single endpoint. If those image pairs that 
include band 4 are deleted, the estimate of single endpoint uncertainty is 

2, 05 pixels. 

It is noted that the pairs that contain band 4 are frequently high. In fact, 
in two instances in scene 3, the pairing was so limited that no lock was 
possible. Loss of lock was also encountered on band 4 for some of the multi- 
temporal measurements (refer to paragraph 3. 3), 

3. 3 Determination of Prefilter and Preprocessor Thresholds 

During the processing of the multispectral image pairs, two peculiarities 
of the LANDSAT image data were observed: 

a, A striping effect between adjacent scan lines, evidently caused by 
uneven sensitivity of the six parallel MSS sensors. The effect is particularly 
noticeable in uniform areas such as water, 

b. A predominance of vertical gradients in the image histograms. 
Although this could be caused in part by striping between scan lines, an 
additional cause might be the tighter spacing of data samples in the horizon- 
tal direction (along scan) than in the vertical direction (between lines). This 
discrepancy, due to the MSS scanning format, would accentuate the vertical 
gradient component. 

The above problems should be relieved by prefiltering of the image data. 
Accordingly, a test was made of the effects of the following filter stages: 
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1. A smoothing filter that averages all square groups of four pixels 
in a uniform manner (single-stage square filter) 

2. Two stages of the above 

3* A filter that averages all vertical groups of two pixels uniformly 
(vertical, correlation factor K = 1) (vertical refers to the crossing of scan 
lines) 

4. A vertical filter with a correlation factor of 0. 4. This value is 
calculated to correspond to the difference between vertical and horizontal 
spacing of pixels in the MSS sensor. 

In addition, it was desired to establish the importance to the measure- 
ments of two preprocessor thresholds, namely; 

• The minimum gradient amplitude that would be accepted 

* The minimum subset length that would be retained. 

The values under consideration for the minimum gradient amplitude were 
2 and 3 gray levels (out of 64), and for the minimum subset length, values 
of 2, 3, and 4 pixels. 

Various combinations of these filters and preprocessor parameters were 
tested on the multitemporal image pair of Figure 2-2, as discussed in 
paragraph 2. 1. The results for relative displacement are shown in Table 
3-2 and plotted in Figure 3-2, 

The columns in Table 3-2 refer to 

a. Type of filter used; One or two sequential stages, using either a 
square (2 x 2) array, or a comparison cietween two pixesl arranged vertically 

b. Minimum detected gradient amplitude 

c. Minimum subset length in pixels 

d. MSS spectral band 

e. The number of subsets obtained 

f. The values for Ax and Ay obtained as described in paragrpah 2. 3, 1 

g. The position of indicated match obtained by linear interpolation 

h. The number of matching endpoints. 
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TABLE 3-2 

EFFECTS OF PREFILTERING AND PREPROCESSOR THRESHOLD SELECTION ON THE 
MEASUREMENTS OF RELATIVE DISPLACEMENT FOR THE WINDOWS IN FIGURE 2-2 


Grad MSL 


Inturpniatton 
X Y 


TwO'Stage 

Square 


TwO'Stagc 

Square 


Two-Stage 

Square 


Two-Stage 

Square 


Two-Stago 

Square 


One-Stage 

Square 


Onc'Slagc 
Vortical 
K = 0.^ 


One-Stage 
Vertical 
K- 1.0 


101/69 

11^/101 

113/1/tO 

85/116 

145/100 
1627152 
161/219 
1 0B/1 65 

138/92 

157/132 

142/204 

113/tSO 


113/112 


197/141 

229/210 


174/129 

203/19D 


273/252 

322/293 


207/179 

221/217 


196/166 

235/222 


0J4-0.39 
0.23 - 0.25 
-0.03 +0.35 
No sign change 

0.07 0.45 
0.26 0.22 
-0.14hQ.13 
-0.14 + 0.47 


0,29-0.25 
0.3 1 -0.20 
No sign change 
No sign change 


0.22 0.43 
0.28-0.10 
-0.33 + 0.17 
-0.25 + 0.12 


0.15-0.34 
0.23 0.14 


0.37 - 0.03 
0.16-0.30 


0.11 -0.36 
0,29 - 0.23 


0.23 - 0.58 
0.30 0.68 
0.24 - 0.24 
0.57 - 0.38 


0.26 - 0.52 
0 38-0.49 
0.41 -O.IQ 
0.65 0.22 

0.35 ■ 0.43 
0.26 0.48 
No sign change 
0.71 0.16 


0.32 0.48 
0.36 - 0.45 
0.22 - 0.19 
0.63 0.44 


0.21 0.62 
0.31 0.36 

0.23 0.46 
134 0.23 


0.24 0.42 
0.29 - 0.35 


0.05 0.44 
0.28 - 0.42 


0 26 

0.28 

■0.48 

0.34 

+0.08 

-0,50 

■0.13 • 

0,33 

0.54 

0.44 

+0:52 

-0.57 

+0.29 

0,75 

-0.54- 

0.45 

■0.61 - 

0.35 

No sign change 

-D.34 

0.40 

-0.74 

0.45 

+0.66 

0.54 

+ 0.68 

0.59 

0.30- 

0.25 

0.62 

0.46 

■0.15 

0.33 

0.0^ 

0.59 

0.51 

0.38 

-0.14 

0.53 

-0.93 

0.36 

0.35 

0.45 

0.23- 

0.10 

-0.56- 

0.40 


7B-0BQ3-T-1 1 








































































Figure 3-2* Effects of Prefiltering and Preprocessor Threshold Selection on 
Relative Displacement Measurements for the Windows in Figure 2-2 








With reference to Figure 3-2, it is apparent that displacements for 
bands 3 and 4 are separated from the results for bands 1 and 2 by an 
average of approximately 0, 5 pixel. However, it is noted from Table 3-2 
that limited numbers of endpoints were available .for these measurements. 
Although ground truth is not available, the cluster of points for bands 1 and 

2 centers at x = -0, 39 and y = -0, 38. All points for these bands are within 
half a pixel of this average position. Since the number of available match 
points is below 100 in half of these measurements, this dispersion is in 
agreement with the discussion of paragraph 3, 2, The effect of the filtering 
and thresholding on the number of available subsets can be seen in Table 

3 " 2, 

The action of the filters in smoothing the gradient histogram is shown 
in Table 3-3. For each of the filter conditions, the number of gradient 
samples for the 16 directions (see Figure 2-4) is tabulated. Both band 1 and 
band 2 are included for both windows of the image pair. For these measure- 
ments, a minimum gradient amplitude of 2 was selected, with a minimum 
subset length of 3, 

At the right of Table 3-3, computations were made of the ratio of the 
number of vertical gradient samples to the number of horizontal samples. 

The vertical filters clearly perform a smoothing action. In a homogeneous 
image, one would expect to achieve values approaching unity for these ratios. 
However, as is evident from Figure 2-5, the horizontal line structure (with 
vertical gradients) predominates over the vertical structure. 

In summary, it was concluded that the choice of a filter and preprocessor 
thresholds .'s not critical, provided that sufficient matching endpoints are 
available. Because it provides smoothing action against sensor striping, as 
well as the calculated amount of vertical correlation, the vertical filter with 
a correlation factor of 0,4 was chosen for the later tests. The selected 
preprocessor thresholds were a minimum gradient amplitude of 2 with a 
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TABLE 3-3 
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GRADIENT DIRECTION HISTOGRAMS FOR VARIOUS PREFILTERS 

FOR WINDOWS IN FIGURE 2-2 



Direction of GraiJItm 
Across Edge 

Ltjihl ^ Dark '^4^ 


1 
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in,j 

lOj) 

13 

9 


Hj 

M 0^^n2 

Grati 

MBL 

Bond 

Impau 
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2 

3 

4 

5 

6 

7 

8 

rr 

10 

n 

12 

13 

14 

IS 

16 

\ 

1. 

No Fitter 

2 

3 

1 

t 

1G3 

204 

2S0 

020 

771 

525 

329 

350 

170 

200 

230 

4B4 

634 

430 

338 

333 

771/169 

631/170 

4.6 

3.7 

8.3 

2. 

One Stiiaa Vert K 0.4 

2 

3 

1 

1 

15^ 

224 

234 

4C7 

530 

303 

217 

4 ID 

171 

210 

216 

330 

447 

330 

300 

412 

5397154 

447/171 

3.5 

2.6 

e.i 

3. 

One-Stage Squari! 

2 

3 

t 

1 

105 

B4 

119 

450 

537 

206 

249 

202 

146 

99 

135 

313 

309 

277 

241 

308 

5347105 

390/146. 

G.1 

2.73 

7.B 

A. 

Tw/o Singe Square 

2 

3 

1 

1 

5B 

70 

122 

432 

400 

227 

232 

244 

104 

72 

101 

279 

337 

236 

210 

253 

4BQ75G 

337/104 

8.7 

3.24 

11.9 

5. 

One-Siage Vml K 1 

2 

3 

1 

1 

HQ 

170 

233 

401 

451 

322 

323 

3D6 

171 

101 

164 

329 

396 

3QD 

254 

420 

451/179 

3057711 

ZS 

2.2 

4.7 

G 




1 

S 
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133 

305 

4D4 

873 

39D 

416 
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1G8 

175 
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389 

628 

399 

204 

207 

813/136 

520/166 

G.3 

3.1 

9.4 

7 




1 

6 
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2ia 

140 

301 

133 

242 

305 

406 
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178 
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307 
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254 

208 
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4337130 

377/153 

3,3 

2.1 

54 

Q 
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1 

5 

7Q 

49 
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307 

420 

210 
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60 

1QG 
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35Q 

2G2 
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236 

429/70 

330/137 

6.1 

2.5 

8,6 

9 




1 

5 

50 

32 

73 
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395 

140 

201 

279 
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bi 

77 
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221 
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2.0 
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mLaimum subset Length of 3. Furthermore, because of Limitations on the 
avaiiabLe number of subsets for this window size, it was decided to perform 
the remaining tests using bands 1 and 2. 

3. 4 Measurement of a Distortion Matrix between MuLtitemporaL Images 
Foliowing the selection of processor parameters, it was decided to per- 
form relative registration measurements on a regular matrix of positions 
between two multitemporal images. The images selected were 1703-17590 
(Reel 2) and 1739-17575 (Reel 2), as shown in Figure 2-2. The matrix con- 
sisted of 16 locations spaced by 200 pixels in both the x and y directions. 

The upper left corner of the first window of the matrix was located on the 
NASA tapes as follows: 

1703-17590: x = 50 

y = 220 

1739-17575: x = 6 

y = 11 

This area corresponds to the mountainous terrain at the top of Figure 2-2. 

In accordance with the discussion introducing paragraph 3. 2, it was de- 
cided to estimate ground truth for these multitemporal images by averaging 
the measurements for bands 1 and 2. The offset will be identical for each 
of these multispectral pairs, but the data for each pair should be largely 
independent. As explained in paragraph 3, 3, use of bands 3 and 4 was re- 
jected because of a limited availability of match points. Although conven- 
tional correlation techniques might be used as an alternate measurement, its 
superior reliability is not assured. Visual comparison between image 
samples could not hope to produce satisfactory results to the precision of 
a small fraction of a pixel as required here. 

Registration measurements for the distortion matrix are contained in 
Tables 3-4 and 3-5. These tables differ in the choice of the match window 
size used for the measurements (see Figure 2-11). For the measurements 
contained in Table 3-4, a match window size of 6 by 6 pixels was used. 
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TABLE 3-4 


MULTITEMPORAL REGISTRATION MEASUREMEN 
IN THE IMAGES IN FIGURE 2-2 (MATCH 


TS FOR A MATRIX OF LOCATIONS 
WINDOW 6x6 PIXELS) 


Locfitian 


Intcrpolaiton 


Motchus 


1/1 2/2 
2 

5/5 1/1 


17/17 18/lB 


21/21 22/22 


25/25 26/26 


29/29 30/30 


33/33 3^i/34 


11 

41/41 42/42 
12 

45/45 46/46 


392/373 

400/400 

374/362 

4GO/4QO 

2QB/255 

400/400 

267/209 

400/400 

400/400 

400/400 

400/358 

400/400 

400/400 

400/400 

400/337 

400/400 

406/400 

400/400 

400/302 

400/400 

40D/3&2 

400/400 

400/29B 

400/400 


■0.07 + 0.63 
0.02 H 0.52 

■0.57 ^ 0.09 
0.05 1- 0.19 

0.10 H 0.37 
•0.23 4 0.35 

■0.26 I 0.36 
^0.13 1*0.36 

■0.21 4- 0,24 
•0.30 t 0.12 

0.20 4 0.31 
0.41 I 0.26 

■0.1 H 0.38 
■0.20 4- 0.25 

0.59 4 0.04 
0,53 4 0.04 

-0.45 i 0*06 
0.3/ 4 0.10 

0.534 0.0B 
-0.20 4- 0.34 

0,24 4 0.22 
0,44 4 0.09 

-0,22 + 0.14 
■0,27 + 0.16 


0,06 

40.53 

0.23 

t 0.49 

0.56 

t 0,05 

0.51 

} 0.00 

0.16 

1 0,61 

0.1.1 

i 0.51 

0.09 

1 0.64 

0.01 

1 0,5>1 

0.30 

1 0.17 

-0,30 

( 0.22 

0.17 

1 0.30 

0.27 

4 0.26 

0,02 

I 0.48 

0,15 

1 0.36 

0.30 

1 0.20 

0.37 

4 0.08 

0.36 

1 0,27 

0,26 

) Q.3B 

0.18 

i 0.22 

0.23 

t 0 23 

0.47 

■ 0.01 

-0.02 

1 0.41 

0.35 

i 0.08 

0.37 

10.14 



2.10 4 6.10 
2.Q4 4 6.32 

3.87 4 5.92 
3.74 + 6,00 

5,13 4 5.21 
5.40 4 5.20 

7.42 + 5.12 
7.22 + 5,07 

0.47 + 3.64 
076 + 3,58 

2.39 + 3,36 

2.61 + 3,51 

4,26 ^ 3.04 
4.44 4 3,29 

5.93 4 2.G0 
5.93 4 3.00 

1 + 0.88 + 1,49 
1 4 0.79 + 1,41 

0.87 + 1,45 
1.38+ 1.50 

2,52 4 0,98 
2.83 4 1,05 

4.61 4 0.81 
4.63 4 0.72 
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TABLE 3-5 


MULTITEMPORAL REGISTRATION MEASUREMENTS FOR A MATRIX OF 
LOCATIONS IN THE IMAGES IN FIGURE 2-2 
(MATCH WINDOW 4x4 PIXELS) 


Location 

Band 

Subsets 

X 

AX 

Y 

AY 

interpolation 
X Y 

Matches 

1 

1/T 

1 

39Z/373 

2 

^0.07 + 0,59 

6 

-0,20 + 0,58 

2,11 

+ 6.26 

248 


2/2 

2 

400/400 

2 

^0.08 + 0,53 

6 

-0.26 + 0.4B 

2,13 

+ 6,36 

282 

2 

5/5 

T 

374/362 

3 

-0,43 + 0.15 

5 

•0.S1 +0.10 

3.74 

+ 5.84 

209 


6/6 

2 

400/400 

3 

0,51+0,16 

5 

•0.61 + 0.07 

3.76 

+ 5.90 

262 

3 

9/9 

1 

288/25B 

5 

•0.13 + 0.51 

5 

-0.12 + 0,67 

5.20 

+ 5,1S 

149 


10/10 

2 

400/4G0 

5 

-0.25 + 0.19 

5 

•0.15 + 0.55 

S.B7 

+ 5.21 

253 

4 

13/13 

1 

207/203 

7 

•0.03+0.91 

5 

‘0,06 + 0,46 

7,07 

+ 5,11 

134 


10/14 

2 

4Q0/40Q 

7 

•0.T3 + 0.54 

5 

-0,07 + 0.57 

7.19 

+ 5.11 

266 

5 

17/ 

1 

400/400 

0 

•0.29 + 0.30 

3 

•0,37+0.15 

0.49 

+ 3.91 

196 


18/ 

2 

400/400 

0 

■0.37 + -0.T3 

3 

•0.39 + 0.13 

0.74 

+ 3,75 

240 

6 

21/ 

1 

400/358 

2 

-0.22+0.29 

3 

-0,21 +0,38 

2.43 

+ 3,36 

180 


22/ 

2 

400/4QD 

2 

-0.42 + 0,24 

3 

-0.3B + 0.29 

2,64 

+ 3.57 

272 

7 

25/ 

T 

400/400 

4 

•0.10 + 0.4T 

3 

‘0,04 + 0.52 

4.20 

+ 3,07 

223 


26/ 

2 

400/4G0 

4 

■0.23 + 0.2B 

3 

-0,03 + 0,39 

4.47 

+ 3.07 

229 

8 

29/ 

1 

400/337 

5 

-0.62 + 0.04 

2 

■0,43 + 0,17 

5,94 

+ 2.72 

188 


30/ 

2 

400/400 

5 

-0,51 + 0.07 

2 

-0,47 + 0,03 

5,88 

+ 2,94 

221 

9 

33/ 

1 

400/400 

‘1 

•0.48+0.16 

1 

-0,17+0,31 

•1 +0.75 +1.35 

209 


34/ 

2 

400/400 

*1 

-0.44 + 0.19 

1 

■0,23+0.23 

-1 +0,69 +1,50 

208 

TO 

37/ 

1 

400/382 

1 

-0.08 + 0.49 

1 

-0.17 + 0.55 1 

1,14 

+ 1.24 

189 


38/ 

2 

400/400 

1 

•0.18 + 0.41 

1 

-0,19+0.27 : 

1.30 

+ 1.41 

22a 

IT 

41/ 

1 

400/392 

3 

•0,00 + 0.47 

0 

-0,64 -0,03 

3,00 

+ 1.00 

175 


42/ 

2 

400/400 

3 

•0.04 + 0.47 

1 

-0.06 + 0.54 

3.0a 

+ 1.10 

207 

12 

45/ 

T 

400/298 

4 

•0,34 + 0.20 

0 

-0,45 + 0.15 

4.63 

+ 0.75 

137 



2 

400/400 

4 

•0,36 + 0,15 

0 

•0.35 + 0.20 

4,71 

+ 0,64 

210 

'"13"" 

49/ 

1 

400/369 

2 

•0.39 + 0.14 

•1 

+0,25+0.30 

•2+0,74 

•1 +0,46 

209 


BO 

2 

400/400 

-2 

•0.35 + 0.17 

*1 

•0,54 + 0,07 

*2 + 0,66*1+0,11 

236 

T4 

53/ 

1 

398/274 

Q 

-0,05 + 0,43 

•1 

•0.21 + 0,51 

0 + 0,10-1 +0.29 

149 


54/ 

2 

400/400 

0 

-0-02 + 0,46 


•0.35+0,20 

0 + 0,04 • 1 + 0,64 

186 

15 

,57/ 

T 

400/329 

1 

•0.45 + 0.14 

'2 

•0.44 + 0,12 

1,74 

• 2 + 0.79 

172 


58/ 

2 

400/400 

T 

•0.43 + 0.01 

2 

•0.32 + 0.14 

2,00 

• 2 + 0,70 

195 

16 

61 

1 

400/342 

3 

•0.29 + 0.01 

•2 

0.29 + 0.16 

3.96 

•2.66 

171 


62 

2 

400/400 

3 

-0.44 + 0.13 

•2 

-0.37 + 0.2s 

3,77 

*2,60 

275 
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For the measurements in Table 3-5, the window size was reduced to 4 x 4 
pixels. Redaction of the window size, if otherwise satisfactory, reduces 
the computational load in the processor for a given image area. Based upon 
the first 12 measurements, the 4 by 4 match window appeared satisfactory, 
hence the matrix of 16 positions was completed with the 4 by 4 window alcne. 

It will be noted in Tables 3-4 and 3-5 that an upper value of 400 subsets 
occurs on many windows. This value represents the subset storage limita- 
tion in the computer simulation program. If this limit is reached in either 
image, the effective registration area between them is reduced at their lower 
end. However, the registration accuracy is based on the number of end- 
points used in the computation; their specific location is not involved. 

A display of the array of measurements for the 16-point matrix is shown 
by Figure 3-3, From these data, it is apparent that a relative orientation 
and scale difference exists between the images. The rotation is approxi- 
mately 0. 24 degree. The scale difference is 0. 9 percent. 

Estimation of registration accuracy using these data will be carried out 
in paragraph 3* 7. 

3. 5 Measurement of a Distortion Matrix between Processed Multitemporal 
Images 

The next tapes received from NASA consisted of Reels No. 2 for images 
1062-15190 and 1080-15192 (Chesapeake Bay region). These tapes were 
understood to be corrected in registration. They each contain a rectangular 
window of 800 pixels in width by 3200 pixels in height. It was desired to 
perform registration measurements across the entire length of this image 
pair. 

A series of 24 registration measurements were made on the pair, includ- 
ing comparison measurements for bands 1 (4) and 2 (5) at 12 locations. 

Four paix'S of measurements were made at a separation of 200 pixels at the 
top of the images; four in the middle and four at the bottom. Data for these 
measurements are contained in Table 3-6 and plotted in Figure 3-4. 
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Figure 3-3- Multitemporal Registration Measurements for a Matrix of 
Locations in the Images in Figure E-3 
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FOR PROCESSED TAPES OF THE 
REGION 
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Figure 3-4, Multitemporal Measurements for Processed Tapes 
of the Chesapeake Bay Region 

It is evident that the distortion between these images is less than 1 pixel. 
However, based upon the clustering of the measurements for top, middle, 
and bottom, a warping to the right exists in the middle of image B relative 
to image A. 

It is of interest that the selection of the measurement matrix did not 
involve prior examination of the image areas. 

3. 6 Quality Measurements on Specially Processed Multitemporal Images 


This series of measurements is a quality determination of an attempt to 
insert a window of image 1739 into image 1703 (see Figure 2-E). The 
coordinates of the insert in image 1703 are shown in Figure 3-5, along with 
the results obtained. The data are contained in Table 3-7, For this test 
only image data for band 2 (5) were supplied. 

The data of Figure 3-5 show that at the top of the insert its position is 
lowered relative to 1703 (y is positive), while at the bottom, it is raised 
(y is negative). In effect, the insert is compressed by nearly 3 pixels in 
the vertical direction. 











TABLE 3-7 

QUALITY MISASUREMENTS ON PROCESSED MULTITEMPORAL IMAGES 


Window 
Number & 
LocatiDn 

Subsets 

1 

50.750 

400/400 

2 

250.750 

400/400 

3 

450,750 

400/400 

4 

650,750 

400/400 

5 

50,1000 

298/349 

6 

250,1000 

400/400 


450.1000 


8 

650,1000 


g 

50,1250 


10 

250,1250 


11 

450,1250 


12 

650,1250 


400/400 


400/400 


291/400 


225/266 


208/237 


320/295 


AX 


-0.41 +0,17 


■O.OB + 0.22 


•0.42+0.16 1 


-0.36+0.10 1 


-0.11 + 0.35 2 


-0.47 + 0.03 1 


-0.43 I 0.08 . 1 


-0.55 + 0.01 1 


-0.17 + 0.10 


-0.23 + 0.29 


-0,08 + 0.39 


•0.04 + 0.25 


Imerpolaton 


X 


Matches 


-0.39 + 0.31 

-1 +DJ2 

1 +0.56 

180 

•0.34^0.14 

-1 + 0.27 

1+0.71 

199 

-0.30 + 0.28 

-1 + 0.72 

1 + 0.52 

191 

•0.42 +0.09 

-1 + 0.78 

1 + 0.82 

223 

-0,01 + 0.44 

0 + 0.24 

2 + 0.02 

142 

■0,54 + -0,01 

-1 +0.92 

1 + 0.9B 

204 

0.30 + 0.25 

-1 +0.84 

1 + 0.55 

239 

-0.21 + 0.26 

0 +0.02 

1 + 0.45 

237 

-0.42 + 0.10 

-1 +0.49 

-2 + 0,81 

159 

-0,02 h 0.35 

0 +0.44 

-1 + 0.05 

135 

-0,53 + 0.10 

0 + 0.17 

-2 + 0.84 

143 

-0.60 + 0.03 

0 + 0.14 

•2 + 0,95 

188 
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3. 7 Estimate of Registration Accuracy 


Determination of registration accuracy is based upon a comparison with 
ground truth. As explained in paragraph 3. 4, an estimate of ground truth is 
made for the multitemporal images by averaging measurements for bands 
1 and 2. Table 3-8 includes a collection of difference components (Ax, Ay) 
for the measurements of bands 1 and 2 from Tables 3-4, 3-5, and 3-6. 

From these components, the magnitude of the radial registration error for 
either band alone is 

A =V (3-U 

This formula assumes that the measured difference between the two dis- 
placement vectors at a given location results from combining two error 
vectors having a random angular relationship. In accordance with equation 
2-3, an estimate for the value of one of them would be the value for A shown 
by equation 3-1. 

In addition to the computed values for A, Table 3-8 indicates the number 
of matching endpoint pairs obtained for both band 1 and band 2 in each case. 

At the bottom of Table 3-8 are shown the mean errors A , for each set of 
measurements and for the group as a whole. The rms value, A , is also 
computed. In accordance with previous comment, is the maximum 

likelihood estimate of the radial registration error. 

A histogram of these error measurements is shown by Figure 3-6, includ- 
ing the location of the average and rms error values. The following comments 
seem appropriate with respect to these results: 

a. None of the estimated errors exceeds 0, 45 pixel 

b. Both the average and rms values of the error are less than 
0, 20 pixel 

c. Two of the three estimates above 0, 30 pixel are based on limited 
numbers of matching endpoints (fewer than 100), 
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TABLE 3-8 

ERROR COMPUTATIONS FROM EARLIER REGISTRATION MEASUREMENTS 



From Table 3 4 

From Table 3 5 

From Table 3*6 

Location 

AX 

AY 

A 

Matches 

AX 

AY 

A 

Matches 

AX 

ay 

A 

Matclies 

1 

0.06 

0,22 

0.16 

290/348 

0,02 

0.10 

0,07 

248/282 

0.08 

0.07 

0.07 

172/228 

2 

0.13 

0.08 

0.11 

280/335 

0,02 

0.06 

0.04 

209/262 

o.oa 

0.16 

0.13 

216/310 

3 

0.27 

0.01 

0.19 

185/319 

0.37 

0.06 

0,27 

149/253 

0.08 

0.03 

0.06 

234/318 

4 

0,22 

0.05 

0.14 

163/323 

0.12 

0 

0.08 

134/266 

0,08 

0.06 

0,07 

136/267 

5 

0.29 

0.06 

0,21 

256/308 

0.25 

0.04 

0.18 

196/240 

0.07 

0.28 

0.20 

202/294 

6 

0.22 

0.15 

0.19 

233/337 

0.21 

0.21 

0.21 

180/272 

0 

0,21 

0.15 

159/251 

7 

0.18 

0.25 

0.22 

312/305 

0.27 

0 

0.19 

228/229 

0.02 

0.22 

0,16 

203/274 

8 

0 

0.40 

0.28 

250/291 

0.06 

0.22 

0.16 

188/221 

0.12 

014 

0.13 

343/373 

9 

0.09 

0.08 

0.09 

267/276 

0.06 

0.15 

0,11 

209/208 

0.36 

0.11 

0.27 

38/115 

10 

0.51 

0.05 

0.36 

256/285 

0.16 

0.17 

0.17 

189/228 

0.05 

0.54 

0,38 

45/109 

11 

0.31 

0.07 

0.22 

254/281 

0.08 

0.10 

0.09 

176/207 

0.59 

0.23 

0.45 

88/200 

12 

0.02 

0.09 

0.07 

183/274 

0.08 

0.11 

0.10 

137/210 

0.17 

0.10 

0.14 

101/163 

13 





0.08 

0.35 

0.25 

209/236 





14 





0.06 

0.35 

0.25 

149/186 





15 





0.26 

0.09 

0.19 

172/195 





16 





0.19 

0.06 

0.14 

171/225 
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0.185 




0.156 




0.184 
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0.201 




0.169 




0.218 
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Figure 3-6. Histogram of Estimated Error for 40 Measurements 

d. The registration measurements were made in accordance with a 
fixed positional matrix. No attempt was made to choose matching areas. 

e. It should be noted that the data from Tables 3-4 and 3-5 are 
correlated to the extent that the same preprocessor outputs were used for 
locations 1-12. However, a different matching window size was used 

(6 by 6 for Table 3-4, 4 by 4 for Table 3-5). 
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Figure 4-1. View of Westinghouse Pattern Recognition Laboratory 
showing the Breadboard Digital Image Processor 

devices, video tape recording equipment, and test instrumentation. A con- 
sole on the bench at right supplies the minicomputer program via punched 
paper tape. A configuration of this kind could be programmed to perform 
registration operations. 

4. 2 Programmable Processor (General-Purpose Computer) 

The output of the preprocessor (subsets) will be fed directly to a buffer 
associated with the programmable processor. The processor must then 
accomplish the following functions; 

1. Binning of subsets in accordance with the registration window 
format 

2. Generation of a search array from the stored subsets of the 
reference window. 

3. Proximity search of subset endpoints from the incoming window 
with respect to the reference window 
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4. Check of matching subset endpoints 

5. Computation of indicated and relative displacements. 

The computation of displacements is a relatively minor portion of the com- 
puter load. Consequently, search operations for up to 10 pixels around a 
given point can easily be accommodated. Initial acquisition should be 
possible between LANDSAT images on this basis. 

The specific determination of pixel throughput rates will depend on a 
detailed examination of processor characteristics and the determination of 
the data format, including desired accuracy, spacing of registration windows, 
and image dimensions. However, it is proposed to select a specific set of 
operating conditions so as to provide an indication of the system capability. 

If we assume an image width of c 0 pixels, of any length, and with registra- 
tion windows spaced 200 pixels apart both across and down, then at a 
1-megapixel-per-second processing rate, a set of four windows will be 
processed every 0. 16 second. If it is further assumed that 200 subsets 
(400 endpoints) will be available, on the average, in each window and that 
150 matches will be obtained per window, then the number of equivalent add 
operations required by the processor will be approximately 150, 000. If the 
processor add time is 1 microsecond, the processor will complete the re- 
quired computations in 0. 15 second, or s’‘ ghtly ahead of the available 0. 16 
second. Accuracy will be in the neighborf'ood of 0. 3 pixel. Initial acquisi- 
tion will be possible as well as continued registration. 

Tradeoffs on the above example can be made with regard to accuracy, 
processing speed, and the size and capability of the general-purpose 
processor. As explained in paragraph 2. 3, accuracies approaching 0. 1 
pixel should be obtained by increasing the level of matches to 800 (a factor 
of 5), The processing load would be increased under these conditions by a 
corresponding amount. 
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5. CONCLUSIONS AND RECOMMENDATIONS 


It is concluded that the primary objectives of this program have been met. 
The Westinghouse registration techniques appear capabLe of performing on 
LANDSAT images, on either a multispectral or multitemporal basis, to 
accuracies of a small fraction of a pixel. They are insensitive to the choice 
of registration area, provided terrain detail is available, which suggests 
application to automated systems. They may be implemented at megapixel- 
per-second rates, using commercial minicomputers in combination with a 
special-purpose preprocessor. In this configuration, they could find appli- 
cation as a rapid means for assessing image geometric quality. 

5. 1 Conclusions 

More detailed conclusions are as follows; 

a. Reasonable agreement was found between the estimated accuracy 
(paragraph 2. 3) and the measured accuracy, both multispectral (paragraph 
3. 2) and multitemporal (paragraph 3. 7). 

b. The highest value of estimated registration error was 0. 45 pixel. 
For a series of 40 measurements, both the mean radial error and rms 
radial error were less than 0. 20 pi::el (paragraph 3. 7). 

c. An indication of measurement quality is given by the number of 
matches involved in the computation (paragraph 2, 3), 

d. Accuracy to 0, 1 pixel could be expected by increasing the 
registration window size. 


5-1 


5. 2 Recommendations 


I 

The following recommendations are made: 

a. For implementation of these registration techniques as a quality- 
control tool, a brief study of the desired system configuration should be 
made, including the available general-purpose processor, image formats, 
and speed and accuracy requirements. Acquisition of a preprocessor should 
be considered, 

b. Additional statistical data collection might be of interest, to in- 
clude; 

1, Effects of cloud cover 

2, Effects of water areas 

3, Additional examination of performance with bands 3 and 4 

4, Measurements to establish accuracy to 0. 1 pixel. 
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